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The opening of a critical-fluctuation induced pseudogap (or precursor pseudogap) in the one- 
particle spectral weight of the half-filled two-dimensional Hubbard model is discussed. This pseu- 
dogap, appearing in our Monte Carlo simulations, may be obtained from many-body techniques 
that use Green functions and vertex corrections that are at the same level of approximation. Self- 
consistent theories of the Eliashberg type (such as the Fluctuation Exchange Approximation) use 
renormalized Green functions and bare vertices in a context where there is no Migdal theorem. They 
do not find the pseudogap, in quantitative and qualitative disagreement with simulations, suggest- 
ing these methods are inadequate for this problem. Differences between precursor pseudogaps and 
strong-coupling pseudogaps are also discussed. 



I. INTRODUCTION 



The two-dimensional Hubbard model is one of the key paradigms of many-body Physics and is extensively studied 
in the context of the cuprate superconductors. While there is now a large consensus about the fact that at half-filling 
(n = 1) the ground state has long-range antiferromagnetic (or spin-density wave) order, the route to this low- 
temperature phase is still a matter of controversy when the system is in the weak to intermediate coupling regime. 
In this regime, we know that the Mermin- Wagner theorem precludes a spin-density-wave phase transition at finite 
temperature but the issue of whether there is, or not, a precursor pseudogap at finite temperature in the single-particle 
spectral weight A(kF,u) is still unresolved. Different many-body approaches give qualitatively different answers to 
this pseudogap question. In particular, the widely used self-consistent Fluctuation Exchange Approximation (FLEX) 
H does not find a pseudogap in the d = 2 repulsive Hubbard model for any filling. A study Q of lattices of up 
to L = 128 found that as the temperature is reduced the quasiparticle peak in A{kp,uj) smears considerably while 
remaining maximum at u) = 0, signaling a deviation from the Fermi liquid behavior but no pseudogap. The same 
qualitative answer is found for attractive models. By contrast, the many-body approach that has given to date the 
best agreement with simulations of both static || and imaginary-time quantities [|| concludes to the existence of a 
precursor single-particle pseudogap in the weak to intermediate coupling regime, for both the attractive and repulsive 
d = 2 Hubbard model, whenever the ground state has long-range order. While we will restrict ourselves to the d = 2 
repulsive model at half-filling, our results will be relevant to the more general question of the pseudogap since small 
changes in filling or changes from repulsive to attractive case [Q Q do not generally necessitate fundamental changes 
in methodology. And the question of many-body methodology is one of our main concerns here. Further comments 
on the regime we do not address here, namely the strong-coupling regime, appear in the concluding paragraphs. 

One may think that numerical results have already resolved the pseudogap issue defined above, but this is not so. 
Early Quantum Monte Carlo (QMC) data analytically continued by the Maximum Entropy method concluded that 
precursors of antifcrromagnctism in A(k, oj) were absent at any non-zero temperature in the weak to intermediate 
coupling regime (U < 8t, U is the Coulomb repulsion term and t the hopping parameter) 0. A subsequent study 
in which a singular value decomposition technique was used instead of Maximum Entropy, concluded to the opening 
of a pseudogap in A(kp,ui) at low temperatures fllpf . Each of the two techniques has limitations. The singular 
value decomposition can achieve a better resolution at low frequencies, but we find that the quality of the spectra 
is influenced by the profile function introduced to limit the range of frequencies. Another difficulty is that it leads 
to negative values of A(\s.,uj). As far as Maximum Entropy is concerned, recent advances (T^j, that we will use here, 
have made this method more reliable than the Classic version applied in Ref . |)| . 

In this paper, we address the issue of the pseudogap in the d = 2, n = 1 Hubbard model at weak to intermediate 
coupling, but it will be clear that the general conclusions are more widely applicable. We present QMC results and 
show that the finite-size behavior obtained for A(kp, u) is correctly reproduced by the method of Ref. [Q. We also 
introduce a slight modification of the latter approach that makes the agreement even more quantitative. This many- 
body approach allows us to extrapolate to infinite size and show that the pseudogap persists even in lattices whose 
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FIG. 1. Formally exact diagrammatic representation of the self-energy in the Hubbard model. The square is the fully 
reducible four-point vertex 



sizes are greater than the antiferromagnetic correlation length £, contrary to the statements made earlier j9|. These 
sizes cannot be reached by QMC when the temperature is too low. We confirm that at low enough temperatures, the 
peak at to — at the Fermi wave vector is replaced by a minimum, corresponding to the opening of a pseudogap Q] 
and by two side peaks that are precursors of the Bogoliubov quasiparticles. In contrast, we find that the A{kp ,uj) 
calculated by FLEX on small lattices are qualitatively different from those of QMC and do not have the correct size 
dependence. Since all many-body techniques involve some type of approximation, their reliability should be gauged 
by their capacity to reproduce, at least qualitatively, the Monte Carlo results in regimes where the latter are free 
from ambiguities. We thus conclude that Eliashberg-type approaches such as FLEX are unreliable in the absence of 
a Migdal theorem and that there is indeed a pseudogap in the weak to intermediate coupling regime at half-filling. 
It is likely, but not yet unambiguously proven, that consistency between the Green functions and vertices used in the 
many-body calculation is crucial to obtain the pseudogap. 



II. MANY-BODY APPROACH 

Many-body techniques of the paramagnon type ]l2] | do lead to a pseudogap but they usually have low-temperature 
problems because they do not satisfy the Mermin- Wagner theorem. No such difficulty arises in the approach of Ref. 
j7j. This method proceeds in two stages. In the zeroth order step, the self-energy is obtained by a Hartree-Fock- 
type factorization of the four-point function with the additional constraint that the factorization is exact when all 
space-time coordinates coincide. jl3| Functional differentiation, as in the Baym-Kadanoff approach [Q, then leads 
to a momentum- and frequency-independent irreducible particle-hole vertex for the spin channel that satisfies j^] 
U sp = U (n^ni) I ((n-f) (nj.)). The irreducible vertex for the charge channel is too complicated to be computed exactly, 
so it is assumed to be constant and its value is found by requiring that the Pauli principle in the form (n^,) = (n a ) 

be satisfied. More specifically, the spin and charge susceptibilities now take the form x^, 1 (q) = Xo(<z) _1 — and 

Xch (l) = Xo(?) _1 + 2" w ^ n ^° com P u ted with the Green function G° that contains the self-energy whose functional 
differentiation gave the vertices. This self-energy is constant, corresponding to the Hartree-Fock-type factorization. 
|fl5f The susceptibilities thus satisfy conservations laws, jlj} the Mermin- Wagner theorem, as well as the Pauli principle 
fflj = (n a ) implicit in the following two sum rules 

Xs p (^ = (("T - n if) =n-2 (n T n x ) (1) 

T 



* ch ^ = (( n T + n lf) - n2 = n + 2 ( n 1 n l) - n2 



where n is the density. We use the notation, q = (q,iq n ) and k = (k,ik n ) with iq n and ik n respectively bosonic and 
fermionic Matsubara frequencies. We work in units where kg = 1, U = 1, lattice spacing and hopping t are unity. 
The above equations, in addition to [|| U sp — U {n^n\) / ((n^) (nj,)), suffice to determine the constant vertices U, 
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and U c h- This Two-Particle Self-Consistent approach will be used throughout this paper, unless we refer to FLEX 
calculations. 

Once the two-particle quantities have been found as above, the next step of the approach of Ref. R|, consists 
in improving the approximation for the single-particle self-energy by starting from an exact expression where the 
high-frequency Hartree-Fock behavior is explicitly factored out. One then substitutes in the exact expression the 
irreducible low frequency vertices U sp and U c h as well as G° (k + q) and Xsp(q), Xch(q) computed above. In the original 
approach || the final formula reads 

(fc) = Un_ a + - 1 Yl P-pXsp (q) + U ch Xch (q)]GUk + q). (2) 

1 

Irreducible vertices, Green functions and susceptibilities appearing on the right-hand side of this expression are all 
at the same level of approximation. They are the same as those used in the calculations of Eq.([l]), hence they are 
consistent in the sense of conserving approximations. The resulting self-energy Y^\k) on the left hand-side though 
is at the next level of approximation so it differs from the self-energy entering the right-hand side. 

There is, however, an ambiguity in obtaining the self-energy formula Eq.(^). Within the assumption that only U sp 
and U c h enter as irreducible particle-hole vertices, the self-energy expression in the transverse spin fluctuation channel 
is different. To resolve this paradox, consider the exact formula for the self-energy represented symbolically by the 
diagram of Fig.l. In this figure, the square is the fully reducible vertex r (q, k — k 1 , k + k' — q) . In all the above formu- 
las, the dependence of T on k + k' — q is neglected since the particle-particle channel is not singular. The longitudinal 
version of the self-energy Eq.(||) takes good care of the singularity of T when its first argument q is near (tt, tt) . The 
transverse version does the same for the dependence on the second argument k — k' , which corresponds to the other 
particle-hole channel. One then expects that averaging the two possibilities gives a better approximation for T since it 
preserves crossing symmetry in the two particle-hole channels. Furthermore, one can verify that the longitudinal spin 
fluctuations in Eq.(g|) contribute an amount U {n^n\) /2 to the consistency condition || ^Tr (E^G ) = U (JifTij.) and 
that each of the two transverse spin components also contribute U (n^ni) jl to ^Tr (E^G ) = U (^t n l) • Hence, aver- 
aging Eq.(||) and the expression in the transverse channel also preserves rotational invariance. In addition, one verifies 
numerically that the exact sum rule |7j — J duj' Im [S CT (k,u/)] /tt = U 2 n- a (1 — n- a ) determining the high-frequency 
behavior is satisfied to a higher degree of accuracy. As a consistency check, one may also verify that ±Tr (£(*)(?(*)) 
differs by only a few percent from^Tr (T,^G°) . We will thus use a self-energy formula that we call "symmetric" 

E« (k) = Un^ + 1 1 ^ U spXs P (q) + U chX ch (<?)] G° (k + q). (3) 

\k) is different from so-called Berk-Schrieffer type expressions j^] that do not satisfy [Jt| the consistency condition 
between one- and two-particle properties, ^Tr (SG) = U (n^ni) . 

In comparing the above self-energy formulas with FLEX, it is important to note that the same renormalized vertices 
and Green function appear in both the conserving susceptibilities and in the self-energy formula Eq.(||). In the latter, 
one of the external vertices is the bare U while the other is dressed (U sp or U c h depending on the type of fluctuation 
exchanged). This means that the fact that Migdal's theorem does not apply here is taken into account. This technique 
is to be contrasted with the FLEX approximation where all the vertices are bare ones, as if there was a Migdal theorem, 
while the dressed Green functions appear in the calculation. The irreducible vertex that is consistent with the dressed 
Green function is frequency and momentum dependent, contrary to the bare vertex appearing in the FLEX self-energy 
expression. In this Eliashberg-type self-consistent approach then, the Green functions are treated at a high level of 
approximation while all the vertices are bare, zeroth order ones. In other words, the basic elements of the perturbation 
theory are treated at extremely different levels of approximation. 



III. MONTE CARLO VS MANY-BODY CALCULATIONS 

Our Monte Carlo results were obtained with the determinantal method Q using typically 10 5 Monte Carlo updates 
per space-time point. The inverse temperature is (3 = 5, the interaction strength is U = 4 and periodic boundary 
conditions on a square lattice are used. Other details about the simulations may be found in the captions. Our 
detailed analysis is for the single-particle spectral weight A (k,oj) at the wave- vector k — (0, tt) but other wave vectors 
will also be shown in the last figure of the paper. The Monte Carlo results are influenced by the statistical uncertainty, 
by the systematic error introduced through imaginary-time discretization, At, and by the finite size, L, of the system. 
The two calculations with At = 1/10 in Fig. 2a show that increasing the number of QMC sweeps (smaller a, defined in 
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FIG. 2. For U = 4, /3 = 5, n = 1, k = (0, n), effect of various other calculational parameters: (a) QMC for L — 6. Thick 
dotted line for At = 1/10 and a = 1.4%. The latter is the average of the error on G(r) normalized by G(r) itself. Calculations 
with the same a but for At = 1/5 and 1/16 are also shown. Thin dashed line is for At — 1/10 but a = 4.2% on G(r). (b) 
Thin solid line, real-frequency calculation using Eq.Q, for an infinite system. Also shown, Maximum Entropy inversion of 
G(t) with same a as in Fig. 3 and a smaller one. 



Fig. 2) leads to a more pronounced pseudogap. The same figure also shows calculations with the same a but different 
values of At (systematic error is of order (At) ). For At ~ 1/10, the decrease in pseudogap depth with decreasing 
At becomes less than the accuracy achievable by the maximum entropy inversion. If the pseudogap persists when 
L — > oo at fixed a and fixed At = 1/10 it should be even more pronounced with a larger number of QMC sweeps 
(smaller a). The size analysis needs to be done however in more detail since increasing the system size L at fixed a 
and At leads to a smaller pseudogap, as shown on the top left panel of Fig. 3a. 

It is customary to analytically continue imaginary time QMC using the Maximum Entropy algorithm To 
provide a faithful comparison with the many-body approaches, we use the imaginary-time formalism for these methods 
and analytically continue them for the same number of imaginary-time points, using precisely the same Maximum 
entropy approach as for QMC. While the round off errors in the many-body approaches are very small, it is preferable 
to artificially set them equal to those in the corresponding QMC simulations to have the same degree of smoothing. 
Many-body results from the symmetric self-energy formula Y,( s \ Eq.(|^), for an infinite system are shown in Fig. 2b. 
The thin solid line is a direct real-frequency calculation in the infinite-size limit. Maximum Entropy inversions of 
the L — y oo value of the many-body G (t) shown on the same figure illustrate that with increasing accuracy the 
real-frequency result is more closely approximated. This confirms that Maximum Entropy simply smooths the results 
when artificially large errors are introduced in the analytical results. ]l6[] For this parameter range, the effects are 
appreciable but do not change qualitatively the results. Even the widths of the peaks are not too badly reproduced 
by Maximum Entropy. The error bars are obtained from the Maximum-Entropy Bayesian probability for different 
regularization parameters a. Jll| They are clearly a lower bound. 

In Fig. 3, we show the spectra obtained for three techniques for system sizes L — 4, 6, 8 and 10. The left-hand panel 
is the QMC data, the middle panel is obtained from E^- 1 Eq.(||) while the last panel is for FLEX. The latter results 
for much larger lattices are not much different from those for the 8x8 system. Since G (k,T) = — j e -p^ +1 A (k, to) , 
the nearly flat (T-independent) portion in G (k,T) of the lower right-hand panel leads, in FLEX, to a maximum in 
A(k,ui) at uj = 0, contrary to the Monte Carlo results. By contrast, as can be seen by comparing the middle and 
left panels, the agreement between Eq.(||) and QMC is very good, except for the height of the peaks. The finite-size 
dependence of the pseudogap for both QMC and Eq.(||) is similar: as the size increases, the depth of the pseudogap 
decreases. Some of the finite-size effects are present in the vertices U sp and U c h- 

Fig. 4a compares three results for the L = 6 system: QMC (thick solid line), and the many-body approach of Ref. 



4 



Monte Carlo Many-Body E (s) FLEX 



L=4 




FIG. 3. Size dependent results for various types of calculations for U = 4, /3 = 5, n = 1, k = (0, n), L — 4,6,8, 10 and 
average relative errors a = 3.4% on G(t). Upper panels show A(\c, uS) extracted from Maximum Entropy on G(r) shown on 
the corresponding lower panels. Each G(r) has 50 points, (a) QMC. (b) Many-body using Eq.Q. (c) FLEX. 

|jj using either the symmetric £^ (Eq.(|^), dotted line) or the longitudinal E" (Eq.(||), thin solid line) self-energy 
formulas. In imaginary time, the agreement between QMC and is striking. The position of the peaks in QMC 
also agrees better with the symmetric version E^ s \ Eq.(|^). 

For the lattice sizes where the Monte Carlo data are qualitatively similar to those of Ref. , and hence uncontrover- 
sial, Fig. 3 has shown that there is a many-body approach that gives good agreement with the simulations. Although 
this many-body approach is not rigorous, especially deep in the pseudogap regime where it is mostly an extrapola- 
tion method these tests suggest that it can give an understanding of finite-size effects in QMC data. There are 
two intrinsic lengths that are relevant, namely £ the antifcrromagnetic correlation length, and £th the single-particle 
thermal de Broglie wavelength defined by vp/T. In simulations, £ may be estimated from the momentum-space width 
of the spin structure factor and £ t h from the Fermi velocity estimated from the maxima of A (k, ui) at different wave 
vectors. For j3 = 5, and L = 10 we have £ ~ 3. At the (tt, 0) point, £th essentially vanishes since we are at the van Hove 
singularity, hence the condition L > t; t h is satisfied. If we had £th > L, one would be effectively probing the finite-size 
zero-temperature quantum regime. When the condition L > ^ t h is satisfied, as is the case here, one has access to the 
finite temperature effects we are looking for. Once agreement on the pseudogap in QMC and the analytical approach 
has been established up to the regime ^ t h < L < £, the analytical approach |7j can be used to reach larger lattice sizes 
(such that £th < £ < L) with relatively modest computer effort. In Fig. 4b we show the spectra obtained by Eq.(|^) 
for L — 6 to 64 and then for L = oo (obtained from numerical integration). We see that the size dependence of the 
pseudogap becomes negligible around L — 32 and that the pseudogap is quite sizable even though it is smaller than 
that in the largest size available in QMC calculations (L = 10) . The size dependence of the pseudogap is qualitatively 
similar when the longitudinal form of the self-energy is used. We thus conclude that the pseudogap exists in the 
thermodynamic limit, contrary to the conclusion of Ref. [pi. The increase in QMC noise with increasing system size 
in the latter work may partly explain the different conclusion. 

The last figure, Fig. 5, shows A (k, ui) obtained by Maximum Entropy inversion of Monte Carlo data (left panel) of 
the many-body approach Eq.(|^) (middle panel) and of FLEX (right panel). Using the symmetry of the lattice and 
particle- hole symmetry, A (k, u>) = A (k+ (tt, tt) , — u>), one can deduce from this figure the results for all wave vectors 
of this 8x8 lattice. The detailed agreement between Monte Carlo and the many-body approach is surprisingly good 
for all wave vectors, even far from the Fermi surface. 
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FIG. 4. U = 4, f3 = 5, n = 1, k = (0,7r), and a = 3.2% in Maximum Entropy, (a) For L = 6, thick solid line for QMC, 
and Many Body using two different self-energy formulas: dashed line for symmetric Eq.Q, and thin solid line for longitudinal 
Eq.(|2|). (b) Size dependent results obtained from symmetric version Eq.Q for L — 6, 8, 10, 16, 32, 64 and infinite size (dashed 
line). The size dependence is monotonous. 



IV. DISCUSSION 



There are two interrelated conclusions to our work. First, detailed analysis of QMC results along with comparisons 
with many-body calculations show that there is a pseudogap in the n = 1, d = 2 Hubbard model, contrary to results 
obtained from previous Monte Carlo simulations M and from self-consistent Eliashberg-type methods such as FLEX. 
Second, we have reinforced the case that the many-body methodology described here is an accurate and simple 
approach for studying the Hubbard model, even as we enter the pseudogap regime. While any self-energy formula 
that takes the form, E oc x (?) C° (k + q) will in general extrapolate correctly to a finite zero-temperature gap Jl7| , 
and hence show a pseudogap as long as \ (?) contains a renormalized classical regime, [Q all other approaches we know 
of suffer from the following defects: they usually predict unphysical phase transitions, they do not satisfy as many 
exact constraints and in addition they do not give the kind of quantitative agreement with simulations that we have 
exhibited in Figs. 3 to 5. Reasons why the mathematical structure of FLEX-type approaches fails to yield a pseudogap 
have been discussed before. The same arguments apply to the pseudogap problem away from half-filling and for 
the attractive Hubbard model as well. Q.||,|l8|,|l9| Since in the Hubbard model there is no Migdal theorem to justify 
the neglect of vertex corrections, it is likely, but unprovcn, that to obtain a pseudogap in FLEX-type approaches, one 
would need to include vertex-correction diagrams that are at the same level of approximation as the renormalized 
Green functions. 

The physical origin of the pseudogap in the 2D Hubbard model has been discussed at great length previously (t]] : The 
precursors of antiferromagnetism in A (kp, to) are preformed Bogoliubov quasiparticles that appear as a consequence 
of the influence of renormalized classical fluctuations in two dimensions. They occur only in low dimension when the 
characteristic spin relaxation rate is smaller than temperature and when £/£th > 1- With perfect nesting (or in the 
attractive Hubbard model) they occur for arbitrarily small U. The ground-state gap value [^o| (and corresponding 
single-particle pseudogap energy scale at finite T) depends on coupling in a BCS-like fashion. 

The previous results show that strong-coupling local particle-hole pairs are not necessary to obtain a pseudogap. 
Such local particle-hole pairs are a different phenomenon. They lead to a single-particle Hubbard gap well above the 
antiferromagnetically ordered state, in any dimension but only when U is large enough, in striking contrast with the 
precursors discussed in the present paper. The Hubbard gap also can exist without long-range order. pl| 

From a methodological point of view, the strong-coupling Hubbard gap is well understood, in particular within 
the dynamical mean-field theory p2] or in strong-coupling perturbation expansion |23[. However, the precursors of 
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FIG. 5. Single particle spectral weight A(k, u>) for U = 4, j3 = 5, n = 1, and all independent wave vectors k of an 8 x 8 
lattice. Results obtained from Maximum Entropy inversion of QMC data on the left panel and many-body calculations with 
Eq.(0) on the middle panel and with FLEX on the right panel. (Relative error in all cases is about 0.3%) 



Bogoliubov quasiparticlcs discussed in the present paper are unobservablc in infinite dimension, where dynamical 
mean- field theory is exact, because they are a low dimensional effect. It remains to be shown if 1/d expansions or 
other extensions of infinite-dimensional methods will succeed in reproducing our results. pl[ 

Experimentally, one can distinguish a strong-coupling pseudogap from a precursor pseudogap (superconducting or 
antiferromagnetic) as follows. Ideally, if one has access experimentally to the critical quantity (spin or pair fluctuations) 
the difference between the two phenomena is clear since precursors occur only in the renormalized classical regime 
of these fluctuations. If one has access only to A (k, lu) , there are also characteristic signatures. The precursors 
are characterized by a "dispersion relation" that is qualitatively similar to that in the ordered state. (However the 
intensity of the peaks in A(k, uj) does not have the full symmetry of the ordered state). By contrast, a strong- 
coupling pseudogap does not show any signs of the symmetry of the ordered state at high enough temperature. |2l| ] 
Also, the temperature dependence of both phenomena is very different since precursors of Bogoliubov quasiparticles 
disappear at sufficiently high t emp erature in a manner that is strongly influenced by the Fermi velocity because of the 
condition £/ (vp/T) > 1. |7|, |l9| , p5f Hence, even with isotropic interactions, the precursor pseudogaps appear at higher 
temperatures on points of the Fermi surface that have smaller Fermi velocity, even in cases when the zero temperature 
value of the gap is isotropic. This has been verified by QMC calculations for the attractive Hubbard model. ^ By 
contrast, at sufficiently strong coupling, the Hubbard gap does not disappear even at relatively large temperatures, 
despite the fact that A (k, u>) may rearrange over frequency ranges much larger than temperature. |2q| 

The methods we have presented here apply with only slight modifications to the attractive Hubbard model case 
where superconducting fluctuations ]l8| may induce a pseudogap [B[l9| in the weak to intermediate coup ling regime 
relevant for the cuprates at that doping |p7f . Recent time-domain transmission spectroscopy experiments [p8[ suggest 
that the renormalized classical regime for the superconducting transition in high-temperature superconductors has 
been observed. Concomitant peaks observed in photoemission experiments p9| persist above the transition tem- 
perature in the normal state. They may be precursors of superconducting Bogoliubov quasiparticles. |8| At exactly 
half-filling on the other hand, the paramagnetic state exhibits a strong-coupling (local particle-hole pairs) Hubbard 
gap. 
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